getwd()
## [1] "/Users/alexg/R files/hair_cortisol/skew-normal FINAL"
library(dlookr)
## Registered S3 methods overwritten by 'dlookr':
## method from
## plot.transform scales
## print.transform scales
##
## Attaching package: 'dlookr'
## The following object is masked from 'package:base':
##
## transform
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(ggplot2)
library(psych)
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
## The following object is masked from 'package:dlookr':
##
## describe
library(reshape)
##
## Attaching package: 'reshape'
## The following object is masked from 'package:dplyr':
##
## rename
library(readxl)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.4 ✔ tidyr 1.3.1
## ✔ readr 2.1.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ psych::%+%() masks ggplot2::%+%()
## ✖ psych::alpha() masks ggplot2::alpha()
## ✖ tidyr::expand() masks reshape::expand()
## ✖ tidyr::extract() masks dlookr::extract()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ reshape::rename() masks dplyr::rename()
## ✖ lubridate::stamp() masks reshape::stamp()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(vioplot)
## Loading required package: sm
## Package 'sm', version 2.2-6.0: type help(sm) for summary information
##
## Attaching package: 'sm'
##
## The following object is masked from 'package:dlookr':
##
## binning
##
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(vtable)
## Loading required package: kableExtra
##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rows
library(brms)
## Loading required package: Rcpp
## Loading 'brms' package (version 2.22.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
##
## Attaching package: 'brms'
##
## The following object is masked from 'package:psych':
##
## cs
##
## The following object is masked from 'package:stats':
##
## ar
library(bayesplot)
## This is bayesplot version 1.12.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
## * Does _not_ affect other ggplot2 plots
## * See ?bayesplot_theme_set for details on theme setting
##
## Attaching package: 'bayesplot'
##
## The following object is masked from 'package:brms':
##
## rhat
library(bayestestR)
library(loo)
## This is loo version 2.8.0
## - Online documentation and vignettes at mc-stan.org/loo
## - As of v2.0.0 loo defaults to 1 core but we recommend using as many as possible. Use the 'cores' argument or set options(mc.cores = NUM_CORES) for an entire session.
library(priorsense)
library(rethinking)
## Loading required package: cmdstanr
## This is cmdstanr version 0.8.0
## - CmdStanR documentation and vignettes: mc-stan.org/cmdstanr
## - CmdStan path: /Users/alexg/.cmdstan/cmdstan-2.36.0
## - CmdStan version: 2.36.0
## Loading required package: posterior
## This is posterior version 1.6.1
##
## Attaching package: 'posterior'
##
## The following object is masked from 'package:bayesplot':
##
## rhat
##
## The following object is masked from 'package:dlookr':
##
## entropy
##
## The following objects are masked from 'package:stats':
##
## mad, sd, var
##
## The following objects are masked from 'package:base':
##
## %in%, match
##
## Loading required package: parallel
## rethinking (Version 2.42)
##
## Attaching package: 'rethinking'
##
## The following object is masked from 'package:loo':
##
## compare
##
## The following objects are masked from 'package:brms':
##
## LOO, stancode, WAIC
##
## The following object is masked from 'package:purrr':
##
## map
##
## The following objects are masked from 'package:psych':
##
## logistic, logit, sim
##
## The following object is masked from 'package:stats':
##
## rstudent
df <- read_xlsx("hair_cort_dog_all.xlsx", col_types = c("text", "text",
"text", "text", "text", "text",
"text", "numeric","text", "skip",
"numeric", "skip", "skip",
"numeric", "skip"))
df <- as.data.frame(df)
dim(df) # will tell you how many rows and columns the dataset has
## [1] 73 11
class(df) # will tell you what data structure has the dataset been assigned
## [1] "data.frame"
head(df)
## number group visit season breed_group coat_colour sex age comorbidity
## 1 c1 stopped v0 winter ret dark Male 43 yes
## 2 c2 stopped v0 autumn mix dark Male 105 yes
## 3 c3 stopped v0 spring ckcs mix Female 117 yes
## 4 c4 stopped v0 summer ret dark Female 108 yes
## 5 c5 stopped v0 summer ret dark Female 110 yes
## 6 c6 stopped v0 winter mix light Female 120 yes
## fat_percent cortisol
## 1 52.21393 4.924220
## 2 38.52059 7.304202
## 3 46.94916 1.590000
## 4 44.46813 0.861570
## 5 39.59363 6.217317
## 6 NA 4.426785
numeric_df <- Filter(is.numeric, df)
describe(numeric_df) # the describe function in psych provides summary stats
## vars n mean sd median trimmed mad min max range skew
## age 1 73 95.85 35.58 101.00 96.07 28.17 16.00 182.00 166.00 -0.10
## fat_percent 2 55 40.48 7.82 40.86 40.70 5.71 17.86 61.20 43.34 -0.28
## cortisol 3 73 8.11 16.47 2.33 4.07 2.17 0.41 104.62 104.20 3.89
## kurtosis se
## age -0.17 4.16
## fat_percent 0.77 1.05
## cortisol 16.80 1.93
plot_normality(numeric_df)
apply(numeric_df, 2, shapiro.test)
## $age
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.97361, p-value = 0.1288
##
##
## $fat_percent
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.97956, p-value = 0.4692
##
##
## $cortisol
##
## Shapiro-Wilk normality test
##
## data: newX[, i]
## W = 0.46269, p-value = 6.756e-15
qqnorm(df$cortisol)
qqline(df$cortisol, col = "red")
qqnorm(log(df$cortisol))
qqline(log(df$cortisol), col = "red")
shapiro.test(log(df$cortisol))
##
## Shapiro-Wilk normality test
##
## data: log(df$cortisol)
## W = 0.94725, p-value = 0.004126
summary(df$cortisol)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4141 1.4119 2.3331 8.1089 6.8455 104.6172
sd(df$cortisol)
## [1] 16.47372
summary(log(df$cortisol))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.8817 0.3449 0.8472 1.1816 1.9236 4.6503
df$lgCort <- log(df$cortisol)
summary(df$lgCort)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.8817 0.3449 0.8472 1.1816 1.9236 4.6503
sd(df$lgCort)
## [1] 1.208074
hist(df$lgCort)
describeBy(cortisol ~ comorbidity, data = df)
##
## Descriptive statistics by group
## comorbidity: no
## vars n mean sd median trimmed mad min max range skew kurtosis
## cortisol 1 15 2.31 1.51 1.84 2.17 0.96 0.75 5.65 4.9 1.19 0.24
## se
## cortisol 0.39
## ------------------------------------------------------------
## comorbidity: yes
## vars n mean sd median trimmed mad min max range skew kurtosis
## cortisol 1 58 9.61 18.2 2.82 5.39 2.89 0.41 104.62 104.2 3.4 12.58
## se
## cortisol 2.39
describeBy(lgCort ~ comorbidity, data = df)
##
## Descriptive statistics by group
## comorbidity: no
## vars n mean sd median trimmed mad min max range skew kurtosis se
## lgCort 1 15 0.66 0.6 0.61 0.65 0.59 -0.29 1.73 2.02 0.23 -0.85 0.16
## ------------------------------------------------------------
## comorbidity: yes
## vars n mean sd median trimmed mad min max range skew kurtosis
## lgCort 1 58 1.32 1.29 1.04 1.24 1.35 -0.88 4.65 5.53 0.59 -0.24
## se
## lgCort 0.17
diffLgCortCo <- mean(df$lgCort[df$comorbidity == "yes"]) - mean(df$lgCort[df$comorbidity == "no"])
diffLgCortCo
## [1] 0.6535362
exp(diffLgCortCo)
## [1] 1.922327
exp(0.25)
## [1] 1.284025
diffLgCortCo / 1.29
## [1] 0.5066172
df$weight_loss_success <- df$group
df$weight_loss_success <- ifelse(df$weight_loss_success == "completed", "yes", "no")
df$breed <- df$breed_group
df$breed <- factor(df$breed, levels = c("mix", "ckcs", "pug", "ret", "other"))
head(df$breed)
## [1] ret mix ckcs ret ret mix
## Levels: mix ckcs pug ret other
df$season <- factor(df$season, levels = c("spring", "summer", "autumn", "winter"))
head(df$season)
## [1] winter autumn spring summer summer winter
## Levels: spring summer autumn winter
df$coat_colour <- factor(df$coat_colour, levels = c("light", "mix", "dark"), ordered = FALSE)
head(df$coat_colour)
## [1] dark dark mix dark dark light
## Levels: light mix dark
sumtable(df)
| Variable | N | Mean | Std. Dev. | Min | Pctl. 25 | Pctl. 75 | Max |
|---|---|---|---|---|---|---|---|
| group | 73 | ||||||
| … completed | 42 | 58% | |||||
| … stopped | 31 | 42% | |||||
| visit | 73 | ||||||
| … v0 | 52 | 71% | |||||
| … v1 | 21 | 29% | |||||
| season | 73 | ||||||
| … spring | 14 | 19% | |||||
| … summer | 22 | 30% | |||||
| … autumn | 21 | 29% | |||||
| … winter | 16 | 22% | |||||
| breed_group | 73 | ||||||
| … ckcs | 7 | 10% | |||||
| … mix | 16 | 22% | |||||
| … other | 26 | 36% | |||||
| … pug | 7 | 10% | |||||
| … ret | 17 | 23% | |||||
| coat_colour | 73 | ||||||
| … light | 27 | 37% | |||||
| … mix | 16 | 22% | |||||
| … dark | 30 | 41% | |||||
| sex | 73 | ||||||
| … Female | 43 | 59% | |||||
| … Male | 30 | 41% | |||||
| age | 73 | 96 | 36 | 16 | 73 | 117 | 182 |
| comorbidity | 73 | ||||||
| … no | 15 | 21% | |||||
| … yes | 58 | 79% | |||||
| fat_percent | 55 | 40 | 7.8 | 18 | 37 | 45 | 61 |
| cortisol | 73 | 8.1 | 16 | 0.41 | 1.4 | 6.8 | 105 |
| lgCort | 73 | 1.2 | 1.2 | -0.88 | 0.34 | 1.9 | 4.7 |
| weight_loss_success | 73 | ||||||
| … no | 31 | 42% | |||||
| … yes | 42 | 58% | |||||
| breed | 73 | ||||||
| … mix | 16 | 22% | |||||
| … ckcs | 7 | 10% | |||||
| … pug | 7 | 10% | |||||
| … ret | 17 | 23% | |||||
| … other | 26 | 36% |
par(mfrow = c(1,1))
vioplot(cortisol ~ sex, col = "firebrick",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ sex, col = "lemonchiffon",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ breed, col = "firebrick",
data = df)
stripchart(lgCort ~ breed, vertical = TRUE, method = "jitter",
col = "steelblue3", data = df, pch = 20)
par(mfrow = c(1,1))
vioplot(lgCort ~ season, col = "lemonchiffon",
data = df)
stripchart(lgCort ~ season, vertical = TRUE, method = "jitter",
col = "steelblue3", data = df, pch = 20)
par(mfrow = c(1,1))
vioplot(cortisol ~ comorbidity, col = "steelblue",
data = df)
par(mfrow = c(1,1))
vioplot(lgCort ~ comorbidity, col = "steelblue",
data = df)
plot(cortisol ~ age, col = "steelblue3", data = df, pch = 20)
plot(lgCort ~ age, col = "steelblue3", data = df, pch = 20)
df$sFat <- standardize(df$fat_percent)
summary(df$sFat)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -2.89165 -0.43568 0.04814 0.00000 0.56366 2.64873 18
df$sAge <- standardize(df$age)
summary(df$sAge)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.2444 -0.6422 0.1448 0.0000 0.5945 2.4215
df$slgCort <- standardize(df$lgC)
summary(df$slgCort)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -1.7079 -0.6925 -0.2768 0.0000 0.6142 2.8713
hist(df$slgCort)
df2 <- na.omit(df)
# plot
ggplot(data = df2, aes(x = comorbidity)) +
geom_boxplot(aes(y = slgCort),
color = "black", width=0.5, lwd = 1) +
geom_point(aes(y = slgCort),
size = 4, color = "purple",
position = position_jitter(width = 0.2)) +
scale_shape(solid = FALSE) +
theme_bw() +
labs(title = "Log hair cortisol by comorbidity") +
labs(y = paste0("Log Hair Cortisol (standardised)")) +
labs(x = paste0("Comorbidity")) +
theme(axis.title.y = element_text(size=16, face = "bold"),
axis.title.x = element_text(size=16, face = "bold"),
axis.text = element_text(size = 14, face = "bold"),
title = element_text(size=18, face="bold"),
plot.title = element_text(hjust = 0.5))
model <- brm(slgCort ~ comorbidity + sAge + breed + sex + (1 | visit), family = skew_normal(), data = df)
default_prior(slgCort ~ comorbidity + sAge + breed + sex + (1 | visit),
family = skew_normal(),
data = df)
## prior class coef group resp dpar nlpar lb ub
## normal(0, 4) alpha
## (flat) b
## (flat) b breedckcs
## (flat) b breedother
## (flat) b breedpug
## (flat) b breedret
## (flat) b comorbidityyes
## (flat) b sAge
## (flat) b sexMale
## student_t(3, -0.3, 2.5) Intercept
## student_t(3, 0, 2.5) sd 0
## student_t(3, 0, 2.5) sd visit 0
## student_t(3, 0, 2.5) sd Intercept visit 0
## student_t(3, 0, 2.5) sigma 0
## source
## default
## default
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## (vectorized)
## default
## default
## (vectorized)
## (vectorized)
## default
NB skip the sigma parameter as this will be determined by the model
# Set individual priors
prior_int <- set_prior("normal(0, 0.5)", class = "Intercept")
prior_b <- set_prior("normal(0, 1)", class = "b")
prior_b_sex <- set_prior("normal(0, 1)", class = "b", coef = "sexMale")
prior_b_co <- set_prior("normal(0.25, 1)", class = "b", coef = "comorbidityyes")
prior_b_sAge <- set_prior("normal(0, 0.5)", class = "b", coef = "sAge")
prior_sd <- set_prior("normal(0, 1)", class = "sd")
prior_alpha <- set_prior("normal(4, 2)", class = "alpha")
# Combine priors into list
priors2 <- c(prior_int, prior_b, prior_b_sex, prior_b_co, prior_b_sAge, prior_sd, prior_alpha)
Increased adapt_delta >0.8 (0.999 here), as had divergent transitions
set.seed(999)
model2 <- brm(bf(slgCort ~ comorbidity + sAge + breed + sex + (1 | visit),
sigma ~ comorbidity),
family = skew_normal(),
prior = priors2,
data = df,
control=list(adapt_delta = 0.99999, stepsize = 0.001, max_treedepth =15),
iter = 8000, warmup = 2000,
cores = 4,
save_pars = save_pars(all =TRUE),
sample_prior = TRUE)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘MacOSX15.5.sdk’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
## Found more than one class "stanfit" in cache; using the first, from namespace 'rethinking'
## Also defined by 'rstan'
summary(model2)
## Family: skew_normal
## Links: mu = identity; sigma = log; alpha = identity
## Formula: slgCort ~ comorbidity + sAge + breed + sex + (1 | visit)
## sigma ~ comorbidity
## Data: df (Number of observations: 73)
## Draws: 4 chains, each with iter = 8000; warmup = 2000; thin = 1;
## total post-warmup draws = 24000
##
## Multilevel Hyperparameters:
## ~visit (Number of levels: 2)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.34 0.35 0.01 1.31 1.00 8851 9509
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -0.59 0.38 -1.35 0.16 1.00 10808
## sigma_Intercept -0.68 0.23 -1.08 -0.18 1.00 13787
## comorbidityyes 0.65 0.24 0.17 1.11 1.00 15043
## sAge -0.14 0.11 -0.34 0.07 1.00 17882
## breedckcs 0.00 0.40 -0.82 0.76 1.00 18541
## breedpug 0.03 0.38 -0.70 0.78 1.00 12245
## breedret -0.10 0.31 -0.71 0.50 1.00 14011
## breedother 0.20 0.31 -0.39 0.81 1.00 12118
## sexMale 0.06 0.22 -0.36 0.49 1.00 19543
## sigma_comorbidityyes 0.77 0.25 0.24 1.23 1.00 15288
## Tail_ESS
## Intercept 13055
## sigma_Intercept 11930
## comorbidityyes 15476
## sAge 15999
## breedckcs 16878
## breedpug 15716
## breedret 16586
## breedother 15436
## sexMale 16482
## sigma_comorbidityyes 13778
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## alpha 3.13 1.30 0.81 6.00 1.00 16416 13614
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(model2)
Looking for hairy caterpillars
mcmc_plot(model2, type = 'rank_overlay')
Usually better than the compatability intervals given in the summary
draws <- as.matrix(model2)
HPDI(draws[,3], 0.97) # 3rd column is draws for comorbidity
## |0.97 0.97|
## 0.1336725 1.1887329
bayes_R2(model2, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles
## Estimate Est.Error Q1.5 Q50 Q98.5
## R2 0.1338352 0.04771195 0.04447219 0.1301455 0.2476114
loo_R2(model2, probs = c(0.015, 0.5, 0.985)) # 0.015, 0.5, 0.985 are the quantiles
## Warning: Some Pareto k diagnostic values are too high. See help('pareto-k-diagnostic') for details.
## Estimate Est.Error Q1.5 Q50 Q98.5
## R2 -0.06298951 0.07628151 -0.2495848 -0.05938339 0.0881158
checks whether actual data is similar to simulated data.
pp_check(model2, ndraws = 100)
par(mfrow = c(1,1))
pp_check(model2, type = "hist", ndraws = 11, binwidth = 0.25) # separate histograms of 11 MCMC draws vs actual data
pp_check(model2, type = "error_hist", ndraws = 11, , binwidth = 0.25) # separate histograms of errors for 11 draws
pp_check(model2, type = "scatter_avg", ndraws = 100) # scatter plot
pp_check(model2, type = "stat_2d") # scatterplot of joint posteriors
## Using all posterior draws for ppc type 'stat_2d' by default.
## Note: in most cases the default test statistic 'mean' is too weak to detect anything of interest.
pp_check(model2, type = "error_scatter_avg")
## Using all posterior draws for ppc type 'error_scatter_avg' by default.
pairs(model2)
loo_model2 <- loo(model2, moment_match = TRUE)
loo_model2
##
## Computed from 24000 by 73 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -102.6 7.0
## p_loo 10.7 2.3
## looic 205.3 14.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
##
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
First, check the sensitivity of the prior and likelihood to power-scaling. Posterior and posteriors resulting from power-scaling.
powerscale_sensitivity(model2,
variable = c("b_Intercept",
"b_sAge",
"b_breedckcs", "b_breedother", "b_breedpug", "b_breedret",
"b_sexMale",
"b_comorbidityyes"))
## Sensitivity based on cjs_dist
## Prior selection: all priors
## Likelihood selection: all data
##
## variable prior likelihood diagnosis
## b_Intercept 0.046 0.051 -
## b_sAge 0.012 0.092 -
## b_breedckcs 0.019 0.072 -
## b_breedother 0.036 0.064 -
## b_breedpug 0.032 0.076 -
## b_breedret 0.019 0.073 -
## b_sexMale 0.014 0.092 -
## b_comorbidityyes 0.021 0.115 -
check_prior(model2, effects = "all")
## Parameter Prior_Quality
## 1 b_Intercept informative
## 2 b_comorbidityyes informative
## 3 b_sAge informative
## 4 b_breedckcs informative
## 5 b_breedpug informative
## 6 b_breedret informative
## 7 b_breedother informative
## 8 b_sexMale informative
## 9 sd_visit__Intercept informative
prior <- prior_draws(model2)
prior %>% glimpse()
## Rows: 24,000
## Columns: 11
## $ Intercept <dbl> 0.49834141, 1.00992523, -0.02432336, -0.19580010, 0.4…
## $ b_comorbidityyes <dbl> 1.39106604, -1.84689682, 0.45871513, -0.87675773, -0.…
## $ b_sAge <dbl> -0.74349270, -0.46312555, -0.01509912, -0.25089330, -…
## $ b_breedckcs <dbl> -0.24004745, -0.94500266, -1.62020698, -0.65152734, -…
## $ b_breedpug <dbl> 0.27830899, -1.73559112, 1.15388048, -0.39574612, 0.3…
## $ b_breedret <dbl> 0.156605962, 0.009103568, -0.149236077, -1.075677793,…
## $ b_breedother <dbl> 0.99324366, 0.48959724, -1.61858092, -1.03557625, 0.1…
## $ b_sexMale <dbl> 0.085446157, 0.241017120, -1.154937454, -0.105595594,…
## $ Intercept_sigma <dbl> -0.48046410, -2.91969964, -1.46334506, 4.45054685, 0.…
## $ alpha <dbl> 3.15414988, 3.38138373, 4.84185186, 3.29769979, 5.266…
## $ sd_visit <dbl> 0.15417643, 0.28194438, 1.82528139, 0.11231738, 0.321…
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_sAge * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
labs(x = "Age (std)",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(-2, 2)) %>%
mutate(c = Intercept + b_comorbidityyes * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "comorbidity",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
# set seed for repeatability
set.seed(666)
# create a new dataframe with 50 draws each from the prior
prior_50 <- prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_comorbidityyes * a) %>%
mutate(comorbidity = ifelse(a == 0, "no", "yes"))
# plot
ggplot(data = prior_50, aes(x = comorbidity)) +
# add a horizontal line at y=0
geom_hline(yintercept=0, linetype="dashed",
color = "grey", size=1) +
# create boxplot
geom_boxplot(aes(y = c),
color = "black", width=0.75, lwd = 1) +
# change the y-axis limits
ylim(-2, 2.5) +
# add points
geom_point(aes(y = c),
size = 4, color = "firebrick3",
position = position_jitter(width = 0.3)) +
scale_shape(solid = FALSE) +
theme_bw() +
labs(title = "Draws from prior for effect of comorbidity on hair cortisol") +
labs(y = paste0("Log Hair Cortisol (standardised)"),
x = paste("Comorbidity")) +
theme(axis.title.y = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 3 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 3 rows containing missing values or values outside the scale range
## (`geom_point()`).
# set seed for repeatability
set.seed(666)
draws <- as.data.frame(draws)
# create a new dataframe with 50 draws each from the prior
draws_50 <- draws %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_comorbidityyes * a) %>%
mutate(comorbidity = ifelse(a == 0, "no", "yes"))
# plot
ggplot(data = draws_50, aes(x = comorbidity)) +
# add a horizontal line at y=0
geom_hline(yintercept=0, linetype="dashed",
color = "grey", size=1) +
# create boxplot
geom_boxplot(aes(y = c),
color = "black", width=0.75, lwd =1) +
# change the y-axis limits
ylim(-2, 2.5) +
# add points
geom_point(aes(y = c),
size = 4, color = "steelblue3",
position = position_jitter(width = 0.3)) +
scale_shape(solid = FALSE) +
theme_bw() +
labs(title = "Draws from posterior for effect of comorbidity on hair cortisol") +
labs(y = paste0("Log Hair Cortisol (standardised)"),
x = paste("Comorbidity")) +
theme(axis.title.y = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_sexMale * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Sex breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedckcs * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "CKCS breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedother * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Other breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedpug * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Pug breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
set.seed(5)
prior %>%
slice_sample(n = 50) %>%
rownames_to_column("draw") %>%
expand_grid(a = c(0, 1)) %>%
mutate(c = Intercept + b_breedpug * a) %>%
ggplot(aes(x = a, y = c)) +
geom_line(aes(group = draw),
color = "firebrick", alpha = .4) +
geom_point(color = "firebrick", size = 2) +
labs(x = "Pug breed",
y = "log(cort) (std)") +
coord_cartesian(ylim = c(-4, 4)) +
theme_bw() +
theme(panel.grid = element_blank())
NB Uses posterior_predict
# use posterior predict to simulate predictions
ppd <- posterior_predict(model2)
par(mfrow = c(2,2))
vioplot(slgCort ~ comorbidity, data = df, main = "Observed", col = "steelblue3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$comorbidity, main = "PPD", col = "firebrick3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$comorbidity, main = "PPD", col = "firebrick3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$comorbidity, main = "PPD", col = "firebrick3")
par(mfrow = c(2,2))
vioplot(slgCort ~ sex, data = df, main = "Observed", col = "steelblue")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$sex, main = "PPD", col = "firebrick3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$sex, main = "PPD", col = "firebrick3")
vioplot(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ df$sex, main = "PPD", col = "firebrick3")
par(mfrow = c(2,2))
stripchart(slgCort ~ breed, vertical = TRUE, method = "jitter",
col = "steelblue3", data = df, pch = 20, main = "Observed")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter",
col = "firebrick3", data = df, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter",
col = "firebrick3", data = df, pch = 20, main = "PPD")
stripchart(ppd[sample(seq(1, dim(ppd)[1]), 1),] ~ breed, vertical = TRUE, method = "jitter",
col = "firebrick3", data = df, pch = 20, main = "PPD")
plot(conditional_effects(model2), ask = FALSE)
ce <- conditional_effects(model2, effects = "comorbidity", method = "posterior_predict")
ce_df <- ce[[1]][c(1, 9:12)]
ggplot(ce_df, aes(x=comorbidity, y=estimate__, group=1)) +
geom_errorbar(width=.1, aes(ymin=lower__, ymax=upper__), colour=c("#F8766D", "#00BFC4"), linewidth = 2) +
geom_point(shape=21, size=10, fill=c("#F8766D", "#00BFC4")) +
theme_bw() +
labs(title = "Conditional effect of comorbidity on hair cortisol") +
labs(y = paste0("Log Hair Cortisol (standardised)")) +
labs(x = paste0("Comorbidity")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey50", size = 12))
mcmc_plot(model2,
variable = c("b_comorbidityyes",
"b_sAge",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sigma_comorbidityyes"))
mcmc_plot(model2,
variable = c("Intercept",
"b_Intercept", "alpha",
"b_sigma_Intercept",
"b_comorbidityyes",
"b_sAge",
"b_sexMale",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sigma_comorbidityyes"))
model2$fit
## Inference for Stan model: anon_model.
## 4 chains, each with iter=8000; warmup=2000; thin=1;
## post-warmup draws per chain=6000, total post-warmup draws=24000.
##
## mean se_mean sd 2.5% 25% 50% 75%
## b_Intercept -0.59 0.00 0.38 -1.35 -0.83 -0.59 -0.35
## b_sigma_Intercept -0.68 0.00 0.23 -1.08 -0.84 -0.69 -0.53
## b_comorbidityyes 0.65 0.00 0.24 0.17 0.49 0.65 0.81
## b_sAge -0.14 0.00 0.11 -0.34 -0.21 -0.14 -0.07
## b_breedckcs 0.00 0.00 0.40 -0.82 -0.26 0.01 0.27
## b_breedpug 0.03 0.00 0.38 -0.70 -0.23 0.03 0.28
## b_breedret -0.10 0.00 0.31 -0.71 -0.30 -0.10 0.11
## b_breedother 0.20 0.00 0.31 -0.39 -0.01 0.19 0.40
## b_sexMale 0.06 0.00 0.22 -0.36 -0.08 0.06 0.21
## b_sigma_comorbidityyes 0.77 0.00 0.25 0.24 0.61 0.79 0.95
## sd_visit__Intercept 0.34 0.00 0.35 0.01 0.09 0.22 0.46
## alpha 3.13 0.01 1.30 0.81 2.27 3.01 3.88
## Intercept 0.00 0.00 0.24 -0.50 -0.13 0.00 0.13
## Intercept_sigma -0.06 0.00 0.10 -0.24 -0.13 -0.06 0.00
## r_visit[v0,Intercept] 0.01 0.00 0.23 -0.49 -0.08 0.00 0.10
## r_visit[v1,Intercept] -0.01 0.00 0.24 -0.53 -0.10 0.00 0.08
## prior_Intercept 0.00 0.00 0.50 -1.00 -0.34 0.00 0.33
## prior_b_comorbidityyes 0.25 0.01 0.99 -1.71 -0.43 0.26 0.92
## prior_b_sAge -0.01 0.00 0.50 -0.98 -0.35 -0.01 0.33
## prior_b_breedckcs 0.01 0.01 1.00 -1.96 -0.66 0.02 0.68
## prior_b_breedpug -0.01 0.01 0.99 -1.97 -0.68 0.00 0.66
## prior_b_breedret 0.01 0.01 1.00 -1.96 -0.66 0.01 0.68
## prior_b_breedother 0.01 0.01 1.01 -1.99 -0.67 0.00 0.69
## prior_b_sexMale 0.01 0.01 1.00 -1.95 -0.67 0.02 0.69
## prior_Intercept_sigma 0.00 0.03 4.33 -7.78 -1.92 0.00 1.90
## prior_alpha 3.99 0.01 2.01 0.10 2.62 3.99 5.36
## prior_sd_visit 0.80 0.00 0.60 0.03 0.32 0.68 1.15
## lprior -10.72 0.01 0.66 -12.51 -10.99 -10.56 -10.27
## lp__ -110.79 0.03 3.02 -117.59 -112.62 -110.45 -108.61
## z_1[1,1] 0.04 0.01 0.77 -1.49 -0.44 0.04 0.52
## z_1[1,2] -0.03 0.01 0.78 -1.60 -0.52 -0.02 0.46
## 97.5% n_eff Rhat
## b_Intercept 0.16 10764 1
## b_sigma_Intercept -0.18 12556 1
## b_comorbidityyes 1.11 14997 1
## b_sAge 0.07 17782 1
## b_breedckcs 0.76 18389 1
## b_breedpug 0.78 12216 1
## b_breedret 0.50 13993 1
## b_breedother 0.81 12112 1
## b_sexMale 0.49 19605 1
## b_sigma_comorbidityyes 1.23 14658 1
## sd_visit__Intercept 1.31 11664 1
## alpha 6.00 15474 1
## Intercept 0.50 9600 1
## Intercept_sigma 0.14 17418 1
## r_visit[v0,Intercept] 0.53 9566 1
## r_visit[v1,Intercept] 0.51 9924 1
## prior_Intercept 0.99 24016 1
## prior_b_comorbidityyes 2.20 23516 1
## prior_b_sAge 0.98 23872 1
## prior_b_breedckcs 1.95 23665 1
## prior_b_breedpug 1.94 24453 1
## prior_b_breedret 1.99 24317 1
## prior_b_breedother 2.00 23671 1
## prior_b_sexMale 1.97 23978 1
## prior_Intercept_sigma 7.78 23398 1
## prior_alpha 7.91 23876 1
## prior_sd_visit 2.25 24697 1
## lprior -9.97 11825 1
## lp__ -105.87 7797 1
## z_1[1,1] 1.63 13136 1
## z_1[1,2] 1.52 13315 1
##
## Samples were drawn using NUTS(diag_e) at Sat Jul 5 23:54:16 2025.
## For each parameter, n_eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor on split chains (at
## convergence, Rhat=1).
mcmc_plot(model2,
variable = c("Intercept", "prior_Intercept",
"b_Intercept"))
mcmc_plot(model2,
variable = c("b_sigma_Intercept", "prior_Intercept_sigma"))
mcmc_plot(model2,
variable = c("b_comorbidityyes", "prior_b_comorbidityyes",
"b_sAge", "prior_b_sAge",
"b_sexMale", "prior_b_sexMale",
"b_breedckcs", "prior_b_breedckcs",
"b_breedpug", "prior_b_breedpug",
"b_breedret", "prior_b_breedret",
"b_breedother", "prior_b_breedother"))
mcmc_plot(model2,
variable = c("b_comorbidityyes", "prior_b_comorbidityyes")) +
theme_classic() +
labs(title = "Prior vs posterior distribution for comorbidity") +
labs(y = "") +
labs(x = paste0("Possible parameter values")) +
scale_y_discrete(labels=c("prior_b_comorbidityyes" = "Prior", "b_comorbidityyes" = "Posterior"),
limits = c("prior_b_comorbidityyes", "b_comorbidityyes")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
mcmc_plot(model2,
variable = c("b_comorbidityyes", "prior_b_comorbidityyes"),
type = "areas") +
theme_classic() +
labs(title = "Prior vs posterior distribution for comorbidity") +
labs(y = "") +
labs(x = paste0("Possible parameter values")) +
scale_y_discrete(labels=c("prior_b_comorbidityyes" = "Prior", "b_comorbidityyes" = "Posterior"),
limits = c("prior_b_comorbidityyes", "b_comorbidityyes")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
mcmc_plot(model2,
variable = c("b_sAge", "prior_b_sAge"))
mcmc_plot(model2,
variable = c("b_sexMale", "prior_b_sexMale"))
mcmc_plot(model2,
variable = c("b_breedckcs", "prior_b_breedckcs"))
mcmc_plot(model2,
variable = c("b_breedpug", "prior_b_breedpug"))
mcmc_plot(model2,
variable = c("b_breedret", "prior_b_breedret"))
mcmc_plot(model2,
variable = c("b_breedother", "prior_b_breedother"))
posterior <- as.matrix(model2)
mcmc_areas(posterior,
pars = c("b_Intercept",
"b_sigma_Intercept",
"alpha",
"b_comorbidityyes",
"b_sAge",
"b_sexMale",
"b_breedckcs",
"b_breedpug",
"b_breedret",
"b_breedother",
"b_sigma_comorbidityyes"),
prob = 0.75) # arbitrary threshold for shading probability mass
posterior <- as.matrix(model2)
mcmc_areas(posterior,
pars = c("b_comorbidityyes",
"b_sigma_comorbidityyes"),
prob = 0.97) + # arbitrary threshold for shading probability mass
theme_classic() +
labs(title = "Posterior distribution for comorbidity parameters",
y = "Density distribution",
x = "Possible parameter values") +
scale_y_discrete(labels=c("b_comorbidityyes" = "Beta", "b_sigma_comorbidityyes" = "Sigma"),
limits=c("b_sigma_comorbidityyes","b_comorbidityyes")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.
# Focus on describing posterior
hdi_range <- hdi(model2, ci = c(0.65, 0.70, 0.80, 0.89, 0.95))
plot(hdi_range, show_intercept = T)
# Focus on describing posterior
hdi_range <- hdi(model2, ci = c(0.65, 0.70, 0.80, 0.89, 0.97), parameters = "comorbidityyes")
plot(hdi_range, show_intercept = T) +
labs(title = "Posterior distribution for comorbidity") +
labs(y = "Density distribution") +
labs(x = "Possible parameter values") +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5),
axis.text.x = element_text(color = "grey50", size = 12),
axis.text.y = element_text(color = "grey8",size = 12))
(when visit = vO, breed = mix, sAge = 0 (mean), sex = Female)
# create new dataframe which contains results of the 10th dog
new_data <- rbind(df2[10,], df2[10,])
# Now change one category to be different
new_data$comorbidity <- c("yes", "no")
# Visualise df to make sure it has worked
new_data
## number group visit season breed_group coat_colour sex age comorbidity
## 12 c12 stopped v0 autumn ckcs dark Male 92 yes
## 121 c12 stopped v0 autumn ckcs dark Male 92 no
## fat_percent cortisol lgCort weight_loss_success breed sFat
## 12 34.01262 1.001929 0.001927641 no ckcs -0.8268174
## 121 34.01262 1.001929 0.001927641 no ckcs -0.8268174
## sAge slgCort
## 12 -0.1081948 -0.9764601
## 121 -0.1081948 -0.9764601
# Now get mean predictions from the draws of the model
pred_means <- posterior_predict(model2, newdata = new_data)
# Compare difference in means for each breed versus mix
difference <- pred_means[,1] - pred_means[,2]
par(mfrow = c(2,2))
# Examine mean of difference for comorbidity
mean(difference)
## [1] 0.6496641
# View histogram of this
dens(difference)
# Create HPDI
HPDI(difference, 0.97)
## |0.97 0.97|
## -2.028449 3.554991
draws <- as.matrix(model2)
mean(draws[,3] >0)
## [1] 0.994375
HPDI(draws[,3], prob=0.97)
## |0.97 0.97|
## 0.1336725 1.1887329
hypothesis(model2, "comorbidityyes >0")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (comorbidityyes) > 0 0.65 0.24 0.25 1.03 176.78
## Post.Prob Star
## 1 0.99 *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(model2, "comorbidityyes <0")
## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (comorbidityyes) < 0 0.65 0.24 0.25 1.03 0.01
## Post.Prob Star
## 1 0.01
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# create new dataframe which contains results of all dogs
new_data1 <- df2
# Now change one category to be different
new_data1$comorbidity <- c("yes")
# create a second new dataframe which contains results of all dogs
new_data2 <- df2
# Now change one category to be different
new_data2$comorbidity <- c("no")
# Now get predictions from the draws of the models
pred_nd1 <- posterior_predict(model2, newdata = new_data1)
pred_nd2 <- posterior_predict(model2, newdata = new_data2)
pred_diff <- pred_nd1 - pred_nd2
pred_diff <- data.frame(pred_diff)
# Create mean of differences for each column (dog) of the dataframe
pred_diff_mean <- apply(pred_diff, 2, mean)
# View histogram of mean differences
hist(pred_diff_mean)
# mean difference
mean(pred_diff_mean)
## [1] 0.6459975
# Create HPDI
HPDI(pred_diff_mean, 0.97)
## |0.97 0.97|
## 0.6350587 0.6582279
pred_slgCort <- posterior_epred(model2)
av_pred_slgCort <- colMeans(pred_slgCort)
plot(av_pred_slgCort ~ df$slgCort)
set.seed(666)
nd <- tibble(comorbidity = c('no', 'yes'), visit = "v0", sAge = 0, sex = "Female", breed = "mix", )
p1 <-
predict(model2,
resp = "slgCort",
newdata = nd) %>%
data.frame() %>%
bind_cols(nd) %>%
ggplot(aes(x = comorbidity, y = Estimate, ymin = Q2.5, ymax = Q97.5)) +
geom_linerange(aes(ymin = Q2.5, ymax = Q97.5),
linewidth = 2, color = c("#F8766D", "#00BFC4"), alpha = 3/5) +
geom_point(size = 8, color = c("#F8766D", "#00BFC4")) +
theme_bw() +
labs(title = "Log hair cortisol (standardised)") +
labs(y = paste0("Log hair cortisol (standardised)")) +
labs(x = paste0("Comorbidity")) +
theme(axis.title.y = element_text(size=12, face="bold"),
axis.title.x = element_text(size=12, face="bold"),
title = element_text(size=12, face="bold"),
plot.title = element_text(hjust = 0.5))
plot(p1)
citation()
## To cite R in publications use:
##
## R Core Team (2025). _R: A Language and Environment for Statistical
## Computing_. R Foundation for Statistical Computing, Vienna, Austria.
## <https://www.R-project.org/>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2025},
## url = {https://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.
citation("dlookr")
## To cite package 'dlookr' in publications use:
##
## Ryu C (2024). _dlookr: Tools for Data Diagnosis, Exploration,
## Transformation_. R package version 0.6.3,
## <https://CRAN.R-project.org/package=dlookr>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {{dlookr}: Tools for Data Diagnosis, Exploration, Transformation},
## author = {Choonghyun Ryu},
## year = {2024},
## note = {R package version 0.6.3},
## url = {https://CRAN.R-project.org/package=dlookr},
## }
packageVersion("dlookr")
## [1] '0.6.3'
citation("dplyr")
## To cite package 'dplyr' in publications use:
##
## Wickham H, François R, Henry L, Müller K, Vaughan D (2023). _dplyr: A
## Grammar of Data Manipulation_. R package version 1.1.4,
## <https://CRAN.R-project.org/package=dplyr>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {dplyr: A Grammar of Data Manipulation},
## author = {Hadley Wickham and Romain François and Lionel Henry and Kirill Müller and Davis Vaughan},
## year = {2023},
## note = {R package version 1.1.4},
## url = {https://CRAN.R-project.org/package=dplyr},
## }
packageVersion("dplyr")
## [1] '1.1.4'
citation("ggplot2")
## To cite ggplot2 in publications, please use
##
## H. Wickham. ggplot2: Elegant Graphics for Data Analysis.
## Springer-Verlag New York, 2016.
##
## A BibTeX entry for LaTeX users is
##
## @Book{,
## author = {Hadley Wickham},
## title = {ggplot2: Elegant Graphics for Data Analysis},
## publisher = {Springer-Verlag New York},
## year = {2016},
## isbn = {978-3-319-24277-4},
## url = {https://ggplot2.tidyverse.org},
## }
packageVersion("ggplot2")
## [1] '3.5.2'
citation("psych")
## To cite package 'psych' in publications use:
##
## William Revelle (2025). _psych: Procedures for Psychological,
## Psychometric, and Personality Research_. Northwestern University,
## Evanston, Illinois. R package version 2.5.3,
## <https://CRAN.R-project.org/package=psych>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {psych: Procedures for Psychological, Psychometric, and Personality Research},
## author = {{William Revelle}},
## organization = {Northwestern University},
## address = {Evanston, Illinois},
## year = {2025},
## note = {R package version 2.5.3},
## url = {https://CRAN.R-project.org/package=psych},
## }
packageVersion("psych")
## [1] '2.5.3'
citation("reshape")
## To cite reshape in publications, please use:
##
## H. Wickham. Reshaping data with the reshape package. Journal of
## Statistical Software, 21(12), 2007.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## author = {Hadley Wickham},
## journal = {Journal of Statistical Software},
## number = {12},
## title = {Reshaping data with the reshape package},
## url = {https://www.jstatsoft.org/v21/i12/},
## volume = {21},
## year = {2007},
## }
packageVersion("reshape")
## [1] '0.8.9'
citation("readxl")
## To cite package 'readxl' in publications use:
##
## Wickham H, Bryan J (2025). _readxl: Read Excel Files_. R package
## version 1.4.5, <https://CRAN.R-project.org/package=readxl>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {readxl: Read Excel Files},
## author = {Hadley Wickham and Jennifer Bryan},
## year = {2025},
## note = {R package version 1.4.5},
## url = {https://CRAN.R-project.org/package=readxl},
## }
packageVersion("readxl")
## [1] '1.4.5'
citation("tidyverse")
## To cite package 'tidyverse' in publications use:
##
## Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
## Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller
## E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V,
## Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to
## the tidyverse." _Journal of Open Source Software_, *4*(43), 1686.
## doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {Welcome to the {tidyverse}},
## author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
## year = {2019},
## journal = {Journal of Open Source Software},
## volume = {4},
## number = {43},
## pages = {1686},
## doi = {10.21105/joss.01686},
## }
packageVersion("tidyverse")
## [1] '2.0.0'
citation("vioplot")
## To cite the enhanced vioplot package in publications use:
##
## Daniel Adler, S. Thomas Kelly, Tom Elliott, and Jordan Adamson
## (2025). vioplot: violin plot. R package version 0.5.1
## https://github.com/TomKellyGenetics/vioplot
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {vioplot: violin plot},
## author = {Daniel Adler and S. Thomas Kelly and Tom Elliott and Jordan Adamson},
## year = {2025},
## note = {R package version 0.5.1},
## url = {https://github.com/TomKellyGenetics/vioplot},
## }
##
## Please also acknowledge the original package: citation("vioplot")
packageVersion("vioplot")
## [1] '0.5.1'
citation("vtable")
## To cite package 'vtable' in publications use:
##
## Huntington-Klein N (2024). _vtable: Variable Table for Variable
## Documentation_. R package version 1.4.8,
## <https://CRAN.R-project.org/package=vtable>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {vtable: Variable Table for Variable Documentation},
## author = {Nick Huntington-Klein},
## year = {2024},
## note = {R package version 1.4.8},
## url = {https://CRAN.R-project.org/package=vtable},
## }
packageVersion("vtable")
## [1] '1.4.8'
citation("brms")
## To cite brms in publications use:
##
## Paul-Christian Bürkner (2017). brms: An R Package for Bayesian
## Multilevel Models Using Stan. Journal of Statistical Software, 80(1),
## 1-28. doi:10.18637/jss.v080.i01
##
## Paul-Christian Bürkner (2018). Advanced Bayesian Multilevel Modeling
## with the R Package brms. The R Journal, 10(1), 395-411.
## doi:10.32614/RJ-2018-017
##
## Paul-Christian Bürkner (2021). Bayesian Item Response Modeling in R
## with brms and Stan. Journal of Statistical Software, 100(5), 1-54.
## doi:10.18637/jss.v100.i05
##
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
packageVersion("brms")
## [1] '2.22.0'
citation("bayesplot")
## To cite the bayesplot R package:
##
## Gabry J, Mahr T (2025). "bayesplot: Plotting for Bayesian Models." R
## package version 1.12.0, <https://mc-stan.org/bayesplot/>.
##
## To cite the bayesplot paper 'Visualization in Bayesian workflow':
##
## Gabry J, Simpson D, Vehtari A, Betancourt M, Gelman A (2019).
## "Visualization in Bayesian workflow." _J. R. Stat. Soc. A_, *182*,
## 389-402. doi:10.1111/rssa.12378 <https://doi.org/10.1111/rssa.12378>.
##
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
packageVersion("bayesplot")
## [1] '1.12.0'
citation("bayestestR")
## To cite bayestestR in publications use:
##
## Makowski, D., Ben-Shachar, M., & Lüdecke, D. (2019). bayestestR:
## Describing Effects and their Uncertainty, Existence and Significance
## within the Bayesian Framework. Journal of Open Source Software,
## 4(40), 1541. doi:10.21105/joss.01541
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {bayestestR: Describing Effects and their Uncertainty, Existence and Significance within the Bayesian Framework.},
## author = {Dominique Makowski and Mattan S. Ben-Shachar and Daniel Lüdecke},
## journal = {Journal of Open Source Software},
## doi = {10.21105/joss.01541},
## year = {2019},
## number = {40},
## volume = {4},
## pages = {1541},
## url = {https://joss.theoj.org/papers/10.21105/joss.01541},
## }
packageVersion("bayestestR")
## [1] '0.16.0'
citation("rethinking")
## To cite package 'rethinking' in publications use:
##
## McElreath R (2024). _rethinking: Statistical Rethinking book
## package_. R package version 2.42, commit
## ac1b3b2cda83f3e14096e2d997a6e30ad109eeee,
## <https://github.com/rmcelreath/rethinking>.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {rethinking: Statistical Rethinking book package},
## author = {Richard McElreath},
## year = {2024},
## note = {R package version 2.42, commit ac1b3b2cda83f3e14096e2d997a6e30ad109eeee},
## url = {https://github.com/rmcelreath/rethinking},
## }
##
## ATTENTION: This citation information has been auto-generated from the
## package DESCRIPTION file and may need manual editing, see
## 'help("citation")'.
packageVersion("rethinking")
## [1] '2.42'
citation("loo")
## To cite the loo R package:
##
## Vehtari A, Gabry J, Magnusson M, Yao Y, Bürkner P, Paananen T, Gelman
## A (2024). "loo: Efficient leave-one-out cross-validation and WAIC for
## Bayesian models." R package version 2.8.0,
## <https://mc-stan.org/loo/>.
##
## To cite the loo paper:
##
## Vehtari A, Gelman A, Gabry J (2017). "Practical Bayesian model
## evaluation using leave-one-out cross-validation and WAIC."
## _Statistics and Computing_, *27*, 1413-1432.
## doi:10.1007/s11222-016-9696-4
## <https://doi.org/10.1007/s11222-016-9696-4>.
##
## To cite the stacking paper:
##
## Yao Y, Vehtari A, Simpson D, Gelman A (2017). "Using stacking to
## average Bayesian predictive distributions." _Bayesian Analysis_.
## doi:10.1214/17-BA1091 <https://doi.org/10.1214/17-BA1091>.
##
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
packageVersion("loo")
## [1] '2.8.0'
citation("priorsense")
## To cite the priorsense package and power-scaling sensitivity:
##
## Kallioinen, N., Paananen, T., Bürkner, P-C., Vehtari, A. (2023).
## Detecting and diagnosing prior and likelihood sensitivity with
## power-scaling. Statistics and Computing. 34(57).
## doi:10.1007/s11222-023-10366-5
##
## To cite Pareto-smoothed importance sampling:
##
## Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, Jonah Gabry
## (2024). Pareto smoothed importance sampling. Journal of Machine
## Learning Research 25(72) https://jmlr.org/papers/v25/19-556.html
##
## To cite importance weighted moment matching:
##
## Topi Paananen, Juho Piironen, Paul-Christian Bürkner, Aki Vehtari
## (2021). Implicitly adaptive importance sampling. Statistics and
## Computing, 31(16), 1-19. doi:10.1007/s11222-020-09982-2
##
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
packageVersion("priorsense")
## [1] '1.1.0'